/*
    This file is part of vlbi2.

    vlbi2 is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    vlbi2 is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with vlbi2.  If not, see <http://www.gnu.org/licenses/>.

#========================================================================
#
# Alan Rogers
# Massachusetts Institute of Technology Haystack Observatory
# Program: vlbi2
# Calculate cross or autocorrelation spectra from Mk5b recorded data.
#
# The vlbi2 program is used to correlate two Mark 5B files with the same
# starting record time and record mode. A postscript file dd1.pos is
# produced showing the correlation amplitude and phase as a function
# of frequency in all channels in the two files with all channels
# plotted on the same axes. Correlation coefficients (fringe amplitudes,
# which will be 1.00 for an autocorrelation) are indicated below each
# channel. The vlbi2 program accepts two files, each 10 MB in size,
# and two flags. The usage is:
#
#        vlbi2 <file1> <file2> [-2bit <2bit_arg>] [-rev <rev_arg>]
#
# Where the items in [] are optional. The -2bit flag indicates whether
# the files are 2 bits/sample recordings (2bit_arg is 1) or 1 bit/sample
# (2bit_arg is 0) [default]. The files must be the same number of
# bits per sample. The .rev flag determines whether the channels in dd1.pos
# are plotted in reverse order (rev_arg is 1) or normal order
# (rev_arg is 0) [default]. It is important to note that with the .rev
# flag, though the channels can be plotted in the reverse order, they
# are not flipped in frequency. If file1 and file2 are the same, an
# autocorrelation is produced; if they are different a cross correlation
# is produced.
#
#
#=======================================================================
#
# Dec 2010
# S.Weston              stuart.weston@aut.ac.nz
# AUT University
#
# Modified to use DISLIN so that we can see the data better than just
# printed across the bottom of an A4 sheet in postscript. Also use
# DISLIN widgets to provide a better UI.
#
#=======================================================================
*/
#include <time.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <stdlib.h>
#include "dislin.h"
#include "vlbi2.h"

int ip, ip1, ip2, id, id_draw, id_print, id_plotallamp,id_plot, id_lis;
int id_plotallphase,id_wgltxt1,id_wgltxt2;
int id_wgltxt3,id_wgltxt4;
int id_amp, id_pha, id_scl;

static char clis[512];
static char *cl1[NCHAN] = {"Channel 1",
                           "Channel 2",
                           "Channel 3",
                           "Channel 4",
                           "Channel 5",
                           "Channel 6",
                           "Channel 7",
                           "Channel 8",
                           "Channel 9",
                           "Channel 10",
                           "Channel 11",
                           "Channel 12",
                           "Channel 13",
                           "Channel 14",
                           "Channel 15",
                           "Channel 16"};

void amp_phase_plot (int id);
void amp_phase_print (int id);
void plot_all_amplitudes();
void plot_all_phases();

void plot(int num,char *file1,char *file2)
{ int i;

  strcpy (clis, cl1[0]);
  for (i = 1; i < NCHAN; i++)
  { strcat (clis, "|");
    strcat (clis, cl1[i]);
  }

  swgtit ("VLBI2 : Correlate Mark5b files");
  swgtyp ("scroll", "list");
  swgpop("NOQUIT");

  swghlp("vlbi2 Version 3 Revision 0 - 14th Dec 2010|This program comes with ABSOLUTELY NO WARRANTY|This is free software, and you are welcome to|redistribute it under certain conditions.|Check gpl-3.0.txt.");

  ip = wgini ("hori");
  swgwth (-20);
  ip1 = wgbas (ip, "vert");
  swgwth (-60);
  ip2 = wgbas (ip, "vert");

/*  swgdrw (2970./2100); */
  swgdrw(2100./2970.); 
  id = wglab (ip1, "Channel:");
  id_lis = wglis (ip1, clis, 1);

  id_plotallamp = wgpbut (ip1, "Plot all amplitudes");
  swgcbk (id_plotallamp, plot_all_amplitudes); 
  id_plotallphase = wgpbut (ip1, "Plot all phases");
  swgcbk (id_plotallphase, plot_all_phases); 

  id_plot = wgpbut (ip1, "Plot selected Channel");
  swgcbk (id_plot, amp_phase_plot); 
  id_print = wgpbut (ip1, "Print");
  swgcbk (id_print,amp_phase_plot);

  id_plot = wgok (ip1);

  id_wgltxt1 = wgltxt (ip1,"File 1:",file1,75);
  id_wgltxt2 = wgltxt (ip1,"File 2:",file2,75);

/* Button widgets to select which plot to produce
   Amplitude
   Phase 
   or both
*/

  id_amp=wgbut(ip1,"Amplitude",1);
  id_pha=wgbut(ip1,"Phase",0);
  id_scl=wgbut(ip1,"Scale equally all channels",0);

  id_wgltxt3 = wgltxt (ip1,"Max Amp","",75);
  id_wgltxt4 = wgltxt (ip1,"Min Amp","",75);

  swgtyp("SCROLL","DRAW");
  id = wglab (ip2, "Phase / Amplitude Plot:");
  id_draw = wgdraw (ip2);

  wgfin ();


}  

void amp_phase_plot (int id)
{

     float amp[NSIZ2],amp_total,amp_avg;
     float y[NSIZ2];
     float phase[NSIZ2];
     int i,ichan,NXA,NYA,NXL,NYL,namp;
     char ctitle[80],c_ampmax[10],c_ampmin[10];
     char cfile[80];
     time_t curtime;

     ichan=gwglis(id_lis);

     if(verbose==1) printf("ichan : %2d\n",ichan);
     if(verbose==1) printf("nsiz2 : %2d\n",NSIZ2);
     ichan=ichan-1;
     amp_total=0.;
     amp_avg=0.;
     namp=0;

     for(i=0;i<NSIZ2;i++) {
        y[i]=i;
/*      Amplitute = sqrt ( sqr real + sqr imaginary ) */

     	amp[i] = sqrt(xspecr[ichan*NSIZ2+i]*xspecr[ichan*NSIZ2+i]
		+xspeca[ichan*NSIZ2+i]*xspeca[ichan*NSIZ2+i]);
      
        amp_avg = amp[i]; 
        amp_total += amp_avg;
        namp++; 

/*      Phase */
	phase[i] = atan2(xspeca[ichan*NSIZ2+i],xspecr[ichan*NSIZ2+i]);

     }

     amp_avg=amp_total/namp;

     if(verbose==1) printf("Average Amplitude : %10.4f\n",amp_avg);

     setxid(id_draw,"widget");

/*
     If this was called by the print button, then setup the metafile for 
     output to file rather than to screen/xwin.
*/
     if(id==id_print) {
         curtime=time(NULL);
         sprintf(cfile,"vlbi2_amplitude_channel%d_%d.pdf",ichan+1,curtime);
         setfil(cfile);
         metafl("pdf");
     } else {

         metafl("xwin");
     }

     setpag("da4p");
     disini();
     erase();
     complx();

     sprintf(c_ampmax,"%f",FindMax(amp));
     swgtxt(id_wgltxt3,c_ampmax);
     sprintf(c_ampmin,"%f",FindMin(amp));
     swgtxt(id_wgltxt4,c_ampmin);

/*   Create Title */
     sprintf(ctitle,"Channel %d",ichan+1);

/*   Check status of button widget for which plot to produce */
     if(gwgbut(id_amp)) {
        if(gwgbut(id_pha)){
          NXA=300;
          NYA=2729;
          NXL=1600;
          NYL=1024;
          axspos(NXA,NYA);
          axslen(NXL,NYL);
        }
        name("Amplitude - total power","Y");
        name("Sample","X");

        if(gwgbut(id_scl)){
             axsscl("LOG","Y");
             setscl(xscl,2,"Y");
        } else {
             setscl(amp,NSIZ2,"Y");
        }

        setscl(y,NSIZ2,"X");
        graf(0.,0.,0.,0.,0.,0.,0.,0.);
        color("red");
        curve(y,amp,NSIZ2);
        color("fore");
        if(gwgbut(id_pha)==0){
           titlin(ctitle,1);
           title();
        }
        endgrf();
     }

     
     if(gwgbut(id_pha)) {
        if(gwgbut(id_amp)){
          NXA=300;
          NYA=1405;
          axspos(NXA,NYA);
        }
        name("Phase","Y");
        name("Sample","X");
        setscl(phase,NSIZ2,"Y");
        axsscl("LIN","Y");
        setscl(y,NSIZ2,"X");
        graf(0.,0.,0.,0.,0.,0.,0.,0.);
        color("green");
        curve(y,phase,NSIZ2);
        color("fore");
        titlin(ctitle,1);
        title();
     }
     disfin();
}

/*=======0========1========2========3========4=======5========6========7========8*/
void plot_all_amplitudes()
{
     float amp[NSIZ2];
     float y[NSIZ2];
     int j,i,ichan,NXA,NYA;
     char message[11];
     char ctitle[80],c_ampmax[10],c_ampmin[10];

     ichan=NCHAN;

/*
     Initialize dislin to plot
*/

     setxid(id_draw,"widget");
     metafl("xwin");
     setpag("da4p");
     disini();
     erase();
     complx();


     NXA=750;
     NYA=50;

     sprintf(c_ampmax,"%f",(float)DFindMax(xspecr));
     swgtxt(id_wgltxt3,c_ampmax);
     sprintf(c_ampmin,"%f",(float)DFindMin(xspecr));
     swgtxt(id_wgltxt4,c_ampmin);

     for(j=0;j<ichan;j++) {

/*      Create Title */
        sprintf(message,"Channel %d",j+1);
 
     	for(i=0;i<NSIZ2;i++) {

       		y[i]=i;

/*      Amplitute = sqrt ( sqr real + sqr imaginary ) */

	        amp[i] = sqrt(xspecr[j*NSIZ2+i]*xspecr[j*NSIZ2+i]
                +xspeca[j*NSIZ2+i]*xspeca[j*NSIZ2+i]);

        }

        axspos(NYA,NXA);
        axslen(500,700);

/*      Was the button to scale equaly set, if yes then plot on log scale
        as values can vary wildly. If the button was not set then scale
        linear
*/
        if(gwgbut(id_scl)){
           axsscl("LOG","Y");
       	   setscl(xscl,2,"Y");
        } else {
       	   setscl(amp,NSIZ2,"Y");
        }

        setscl(y,NSIZ2,"X");
        setgrf("line","line","line","line");
	graf(0.,0.,0.,0.,0.,0.,0.,0.);
        messag(message,NYA+90,NXA-50);
	color("red");
	curve(y,amp,NSIZ2);
	color("fore");
        endgrf();

        NYA=NYA+500;
        if(NYA>1550){
            NYA=50;
            NXA=NXA+700;
        }
     }

     disfin();
}

/*=======0========1========2========3========4=======5========6========7========8*/
void plot_all_phases()
{
     float phase[NSIZ2];
     float y[NSIZ2];
     int j,i,ichan,NXA,NYA;
     char message[11];

     ichan=NCHAN;

/*
     Initialize dislin to plot
*/

     setxid(id_draw,"widget");
     metafl("xwin");
     setpag("da4p");
     disini();
     erase();
     complx();


     NXA=750;
     NYA=50;

     for(j=0;j<ichan;j++) {

/*      Create Title */
        sprintf(message,"Channel %d",j+1);

        for(i=0;i<NSIZ2;i++) {

                y[i]=i;

/*              Phase */
                phase[i] = atan2(xspeca[j*NSIZ2+i],xspecr[j*NSIZ2+i]);

        }

        axspos(NYA,NXA);
        axslen(500,700);
        setscl(phase,NSIZ2,"Y");
        setscl(y,NSIZ2,"X");
        setgrf("line","line","line","line");
        graf(0.,0.,0.,0.,0.,0.,0.,0.);
        messag(message,NYA+90,NXA-50);
        color("red");
        curve(y,phase,NSIZ2);
        color("fore");
        endgrf();

        NYA=NYA+500;
        if(NYA>1550){
            NYA=50;
            NXA=NXA+700;
        }
     }

     disfin();

}

/*=======0========1========2========3========4=======5========6========7========8*/
void amp_phase_print (int id)
{
     printf("***WARNING: Printing not enabled yet\n");
}
